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We describe the construction of stepped-pressure equilibria as extrema of a multi-region, relaxed 
magnetohydrodynamic (MHD) energy functional that combines elements of ideal MHD and Taylor 
relaxation, and which we call MRXMHD. The model is compatible with Hamiltonian chaos theory 
and allows the three-dimensional MHD equilibrium problem to be formulated in a well-posed manner 
suitable for computation. The energy-functional is discretized using a mixed finite-element, Fourier 
representation for the magnetic vector potential and the equilibrium geometry; and numerical solu- 
tions are constructed using the stepped-pressure equilibrium code, SPEC. Convergence studies with 
respect to radial and Fourier resolution are presented. 



I. INTRODUCTION 

Zcro-Larmor-radius, single- fluid magnctohydrodynam- 
ics (MHD) is commonly used for modeling the global, 
long-time-scale state of plasmas in the magnetic confine- 
ment devices used for fusion power research. It is often 
reasonable to approximate the plasma pressure tensor as 
isotropic and to ignore inertial effects due to small mass 
flows. There is no minimum length scale in this model, 
so spatial discontinuities are allowed [1]. To allow a weak 
formulation, we write the equilibrium condition in con- 
servation form, 

v-jW^i-^Uo, (i) 

V 2^ Mo / 

where, using SI units, fig is the permeability of free 
space, p(r) > is the pressure as a function of position 
r = xi + yj + zk, and B(r) is the magnetic field, which 
must obey V-B = 0. While MHD is a rather crude model 
for the physics of a plasma, the Maxwell equations for the 
magnetic field and the 'dynamics' of field lines are exact. 

The problem addressed in this paper is, treating both p 
and B as unknown fields within suitable function spaces, 
find general, weak solutions of Eq. (1) in an arbitrary, 
three-dimensional (3D) toroidal domain, V, under the 
homogeneous boundary conditions, 

p = 0, n-B = 0, VredV, (2) 

where n is the unit normal at the boundary, dV. We 
take the boundary to be fixed, being either the edge of a 
plasma confined by a notional, tight-fitting shell, or the 
boundary of a surrounding vacuum region. 

Our goal is to formulate the 3D equilibrium problem 
in a way that is as well-posed mathematically as the two- 
dimensional (2D) problem, by which we mean that a well- 



*Elcctronic address: shudson@pppl.gov 
t Electronic address: robert.dewar@anu.edu.au 
+ Electronic address: graham.dennisOanu.edu.au 
§ Electronic address: matthew.hole@anu.edu.au 
^Electronic address: mathew.mcgannOanu.edu.au 
"Electronic address: greg.vonnessi@anu.edu.au 
^Electronic address: slazersoOpppl.gov 



defined, unique solution exists. And, to develop an ac- 
curate, robust and efficient numerical solution method, 
where the error between the approximate numerical so- 
lution, e.g. fh, and the exact solution, /, is bounded and 
goes to zero as fh = f + 0(h n ), where h characterizes 
the numerical resolution and n depends on the numerical 
discretization. 

If p and B are assumed to be diffcrentiable within a 
subregion of V, then Eq. (1) is locally equivalent to the 
force-balance condition 

Vp=jxB, (3) 

where j = V x B is the current density (here, and here- 
after, fj,Q is ignored). We will not restrict attention to 
diffcrentiable solutions in the following, but we will work 
within the approximation that B ■ Vp = 0, which follows 
directly from Eq. (3). Physically, this approximates the 
transport of heat and mass along the magnetic field as in- 
finite compared to that across the field. This immediately 
implies that p is invariant along the magnetic field: the 
spatial dependence of the pressure and the phase space 
structure of the magnetic field are intimately connected. 

A specific equilibrium state is characterized by the 
pressure, i.e. p is considered to be a supplied, input func- 
tion. The computational challenge is to then determine 
the magnetic field that is consistent with the given pres- 
sure and boundary. Gencrically, in 3D, there exist regions 
within V where the magnetic field lines are chaotic. To 
admit numerically tractable solutions for B, it is neces- 
sary to restrict the class of admissible functions for p: 
and to guarantee that B is consistent with a given p, 
topological constraints on B must be enforced. 

In Sec. II, we review the salient properties of 3D mag- 
netic fields, which generally have a fractal phase space, 
and we sketch the nature of continuous solutions for p and 
B. This is based on the construction of an crgodic par- 
tition; which, being fractal, is impractical from a stand- 
point of numerical implementation. So, we describe a 
discrete partition which greatly simplifies the equilibrium 
problem and leads naturally to stepped-pressure equi- 
libria, where the plasma is modeled as a set of nested 
volumes in each of which the field satisfies the Beltrami 
equation, V x B = ^B, and across the interfaces that 
separate these volumes the total pressure is continuous, 
[b + J B 2 /2]]=0. 
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'Sharp-boundary' [2] states and multi- volume [3] 
sharp-boundary states have been considered previously, 
and Bruno & Laurence [4] have presented theorems that 
insure the existence of sharp boundary solutions, with 
an arbitrary number of pressure jumps, for tori whose 
departure from axisymmetry is sufficiently small. In 
Sec. Ill we introduce a variational approach to solv- 
ing Eq. (1) based on the notion of multi- region, relaxed 
MHD (MRXMHD), which is a generalization of Taylor's 
[5] relaxed-MHD formulation: that a sufficiently tur- 
bulent/chaotic, weakly non-ideal plasma will evolve so 
as to minimize the energy subject to the constraint of 
conserved magnetic helicity, and in doing so will break 
most of the constraints [6] of ideal MHD, thus allow- 
ing magnetic reconnection. In MRXMHD, a plasma 
with a non-trivial pressure profile is constructed as a 
nested collection of relaxed states, between which the 
ideal-MHD constraints apply. By deriving the Euler- 
Lagrange equations, we see that the MRXMHD energy 
functional has stepped-pressure equilibria as extremiz- 
ing solutions. A close examination of the force-balance 
condition, \[p + B 2 /2]] = 0, reveals that the rotational 
transform of the interfaces must be strongly irrational. 

In Sec. IV, the MRXMHD energy functional is dis- 
cretized using a mixed Fourier, finite-clement represen- 
tation for the vector potential and geometry. Setting to 
zero the derivatives of the energy functional with respect 
to the vector-potential in each volume gives a linear sys- 
tem for the magnetic field, V x B = fiR, where /i is a 
Lagrange multiplier (sometimes called the Beltrami pa- 
rameter). This can be adjusted in order to preserve the 
helicity integral, or both /i and the enclosed poloidal flux 
can be adjusted to satisfy the interface rotational trans- 
form constraints. 

Assuming the Beltrami fields in each volume have been 
computed for a arbitrary interface geometry, the prob- 
lem of constructing an equilibrium solution is standard: 
changes in the interface geometry are allowed to cither 
minimize the energy functional using conjugate gradient 
methods, or to find a zero of the multi-dimensional gra- 
dient = force-balance vector using a Newton method. To 
fully constrain the Fourier representation of the interface 
geometry, we employ spectral-condensation [7, 8] meth- 
ods to obtain a preferred poloidal angle coordinate. Il- 
lustration of equilibrium states and convergence studies 
are then presented. 

At appropriate points in the discourse, we contrast our 
approach to constructing equilibrium solutions with oth- 
ers in the literature. 



II. HAMILTONIAN CHAOS, PARTITIONED 

Magnetic- field-line flow is a Hamiltonian system [9]. 
The well-developed theory of Hamiltonian dynamical sys- 
tems (see, for example, the texts by Wiggins [10] and 
Lichtenberg & Lieberman [11], and the review by Mciss 
[12]) provides a strong foundation on which to build. We 
shall sometimes use general dynamical-systems language 
rather than the more specialized plasma terminology; 
for instance, using 'orbit' and 'magnetic field line' in- 



terchangeably. To facilitate the following discussion, we 
use cylindrical coordinates, (R, <$>, Z), which are orthogo- 
nal and right handed, so that x = Rcos(<fi), y = Rsm(<fi), 
and z = Z, and (x, y, z) are Cartesian. 

Devices of the tokamak and reversed-field-pinch (RFP) 
[9] classes use a large number of identical toroidal field 
coils arranged with a discrete rotational symmetry about 
the z-axis. In the axisymmetric special case, it is reason- 
able to seek solutions that are invariant under rotation. 

Axisymmetric magnetic fields are representable as 1- 
degree-of- freedom (1-dof) autonomous Hamiltonian sys- 
tems [9], with <p, periodic, playing the role of time. Such 
systems are integrable in the dynamical systems sense, 
and action-angle coordinates may be constructed. The 
field lines lie on nested invariant tori, %p — const, that 
foliate the extended phase space, (ip,9,(j)), where tp is a 
toroidal flux function and 9 is a poloidal angle that in- 
creases linearly against </>. In the terminology of magnetic 
confinement, the invariant tori are called magnetic flux 
surfaces, and action-angle coordinates are called straight- 
ficld-linc coordinates. In the following, when we refer to 
an integrable system, we will assume that the integrable 
system has shear. 

The invariant tori = flux surfaces are characterized 
by their rotation number = rotational transform, which 
is commonly denoted in magnetic confinement plasma 
physics by t. (Historically [13], the term rotational trans- 
form refers to t, the average poloidal angle increase in 
each iteration of the return map, but in modern usage 
[9] it is used for i = l/2tt, and often the 'bar' is omit- 
ted.) If « is a rational number, i = n/m where n and 
m are integers, then the corresponding surface is foliated 
by periodic orbits = closed field lines, which close on 
themselves after m toroidal transits, having undergone 
n poloidal transits. If t is an irrational number, then 
the flux surface is covered ergodically by a single quasi- 
periodic orbit, which never closes on itself (but comes ar- 
bitrarily close) , and each irrational surface is the closure 
of an irrational field line. 

In the axisymmetric case, the equilibrium problem can 
be reduced to the task of solving a 2D partial differential 
equation, the Grad-Shafranov equation [9, 14], which is 
well-posed (except for bifurcations [15]). The equilibria 
are characterized by two free profile functions, e.g. the 
pressure, p(ip), and the rotational transform, t(ip). Be- 
cause space is foliated by flux surfaces, equilibria with 
continuous, smooth profiles are admissible; in fact, the 
only magnetic fields consistent with B • Vp = and glob- 
ally smooth profiles are integrable magnetic fields. 

Axisymmetry is necessarily always broken to some ex- 
tent by the modular nature of the conductors and ma- 
chine imperfection, or by intentionally applied perturba- 
tion fields [16], or by equilibrium bifurcations [17]. The 
stellarator family [9, 18] of confinement devices is inten- 
tionally nonaxisymmetric. This allows greater freedom 
in the design of experiments and can provide enhanced 
plasma stability. (The nonaxisymmctry of stellarators, 
however, generally leads to degraded particle confine- 
ment; this can be ameliorated, somewhat, by the use of 
'quasi-symmetric' configurations [19].) 

The 3D magnetic field-line flow is still analogous to a 
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Hamiltonian dynamical system, but because there is no 
longer a symmetry coordinate, the 3D field-line Hamil- 
tonian is not autonomous. Such systems, still periodic 
in (j>, are called li-dof systems and are generically non- 
integrable, meaning that the extended phase space is al- 
most never foliated by invariant tori. 

The periodic orbits are fragile. Resonant magnetic 
fields associated with geometric deformation destroy al- 
most all of the periodic orbits; magnetic islands form, 
and regions of chaotic magnetic field lines emerge. The 
destruction of rational surfaces is related to the classical 
problem of small denominators in the transformation to 
action-angle coordinates for the perturbed system. The 
Poincarc-Birkhoff theorem [12] shows that, for every ra- 
tional invariant torus present in the integrable case, at 
least two of the periodic orbits survive. One orbit is hy- 
pcrbolically unstable, while the other is elliptically stable 
or has become hyperbolic through a period-doubling bi- 
furcation. 

These orbits, known as Poincarc-Birkhoff orbits, form 
a robust 'skeleton' of invariant sets and provide crucial 
information about the structure of phase space. The ex- 
istence of a given KAM surface (described below) can be 
inferred from the stability of nearby periodic orbits using 
Greene's residue criterion [20]. Associated with each un- 
stable periodic orbit is an unstable manifold and a chaotic 
sea [10], comprised of irregular trajectories without a well 
defined rotational transform, i.e. the ratio AO/ A<j> does 
not converge as A<p — ¥ oo, where A9 and A</> are the 
increase in 9 and <f> along a field line. Although there is 
no formal proof [12], it is standard to assume, based on 
computational evidence, that the closure of each chaotic 
sea is a three-dimensional subset of R 3 , as each irregular 
trajectory seems to fill a volume. Associated with the 
elliptic periodic orbits are local regions of regular trajec- 
tories, the so-called magnetic islands. 

The irrational field lines are quite robust to pertur- 
bation. Indeed, they are guaranteed to survive by the 
Aubry-Mather theorem; however, a given irrational field 
line may or may not come arbitrarily close to every point 
on a smooth surface. The KAM theorem, named in honor 
of Kolmogorov, Arnold and Moser [21-25], shows that 
a finite measure of invariant tori do exist for sufficiently 
small, smooth perturbations to an integrable system, pro- 
vided that the rotational transform, t, is sufficiently irra- 
tional, i.e. * must satisfy a Diophantine condition: there 
exists an r > and k > 2 such that, for all integers n and 
m, \t — n/m\ > r/m k . About each rational, n/m, there 
is an excluded region of width r/m k , which is consistent 
with the emergence of a chaotic sea about every unstable 
periodic orbit. KAM tori are two-dimensional subsets of 
R 3 whose union is of finite measure and forms a partition 
of phase space. 

Typically, as the magnitude of the geometric deforma- 
tion increases, the size of the magnetic islands increases, 
the volume of the chaotic seas increases, and each given 
KAM surface will become more geometrically deformed 
until a critical point is reached, at which point the sur- 
face is continuous but no longer smooth. These critical 
tori form fractal boundaries between the chaotic seas as- 
sociated with different island chains. By 'fractal' we sim- 



ply mean having a hierarchy of qualitatively self-similar 
structure on all scales, with no minimum length scale, 
and is non-diffcrcntiable. 

Some KAM tori are more robust than others. The 
most robust invariant tori are those that have the 'most 
irrational' rotational transforms, where 'most irrational' 
means most difficult to approximate with rationals. Such 
irrationals arc called noble, and their definition is made 
precise using the continued fraction representation [26]. 
The noble KAM tori are also the smoothest, in that fewer 
Fourier harmonics are required for an accurate descrip- 
tion of their geometry. 

After the destruction of a KAM surface, the closure of 
an irrational field line has the structure of a Cantor set 
[27, 28] and is called a cantorus [29] (hint: Cantor + torus 
= cantorus). Cantori are one-dimensional subsets of R 3 
[30], and constitute a set of zero measure that does not 
serve to partition phase space. The cantori can, however, 
form effective partial barriers to field-line transport [31] 
and thus also to anisotropic diffusion [32] . 

Ergodic invariant sets form a fractal hierarchy. The 
'primary' chaotic seas and KAM tori are infinitely inter- 
twined, and each chaotic sea contains 'secondary' island 
chains, KAM tori, cantori, and chaotic seas in an 'islands 
around islands' pattern repeated ad infinitum [12, 33- 
36]. The chaotic seas are infinitely multiply-connected, 
bounded externally by critical, primary invariant tori, 
and internally by the infinite hierarchy of islands. 



A. continuous solution on ergodic partition 

To understand the general class of functions for p and 
B that admit solutions to the equilibrium problem we 
first consider the implications that the non-integrability 
of the magnetic field has on the structure of the pressure; 
and then, given a pressure that is consistent with a non- 
intcgrablc field, consider the implications this has on the 
field itself. 

To understand the structure of the pressure function 
that satisfies B • Vp = 0, given a generic magnetic field, 
it is convenient to represent phase space as a collection of 
pair-wise disjoint sets that are invariant under the field- 
line flow map, ip<p : ro 4 r. This map is constructed 
simply by following a field line a distance <f> in toroidal 
angle from a point, ro, on a surface of section (for exam- 
ple the <f> = plane) to arrive at point r. The return map 
is generated by following field lines once around the ma- 
chine back to the initial surface of section. (In the case of 
the RFP device, this discussion applies in a subdomain 
not containing points where reverses sign - to treat 
the toroidal field-reversal region a poloidal surface of sec- 
tion should be used instead.) An invariant set icl 2 
within a surface of section is a set invariant under the re- 
turn map, <P27r(^4) = -4- An invariant set V C R 3 within 
phase space may be constructed as the continuous union 

Of SUCh SCtS, i.C. V = U0e[O,27r) V4>( A )- 

An invariant partition is a union of invariant tori, 
which arc two-dimensional magnetic surfaces, and 
three-dimensional invariant toroidal volumes or toroids 
bounded by invariant tori. If an invariant volume con- 
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tains an invariant surface, e.g. a KAM surface, then the 
volume may be subdivided into two distinct subvolumes, 
each of which is invariant under the return map. An er- 
godic invariant set is a set with finite measure that allows 
no further subdivision, and the ergodic partition of phase 
space is its decomposition into ergodic invariant sets and 
a nonergodic (periodic) set of zero measure - see defini- 
tion 2.1 of Ref.[35]. 

For the purpose of constructing weak solutions to the 
equilibrium problem, we are primarily interested in the 
sets of finite measure. We ignore the cantori and periodic 
orbits and take our partition as having every chaotic sea, 
C a , each of which has finite volume, and the invariant 
surfaces, Cp, the union of which has finite measure, as 
the only elements with non-trivial measure, where a and 
j3 are elements of appropriate indexing sets (e.g. a is a 
rational and j3 is irrational) . As any field line approaches 
every point arbitrarily closely in a given ergodic set, the 
only solution for p consistent with B • Vp = is p = p a = 
const, Vr G C a , and similarly for Cp. The most general 
solution for the pressure is 

p(r)=p a l a (r)+pplp(r), (4) 

where I a is an indicator function on each ergodic com- 
ponent, i.e. I a (r) = 1 if r G C a and 7 Q (r) = otherwise, 
and similarly for I p. 

We now recognize that B(r) is not arbitrary and seek 
a similarly general characterization of the constraints 
that Eq. (3) places on this function. Each C a has finite 
volume, and we assume that B is diffcrcntiable within 
C a . That the pressure is constant in C a implies that 
Vp = 0. Then, force balance, Vp = j x B, implies that 
V x B = ^(r)B, for some scalar function /z(r). Taking 
the divergence of this equation, we find B \7fi = 0. Thus, 
like p, \x must be constant within each ergodic region [37] , 
fx = Li a in C a , and B must be a linear force- free field, i.e. 
it satisfies the Beltrami equation, 

VxB = /i Q B. (5) 

This is a well-studied linear elliptic partial differen- 
tial equation, about which much is known [38-43]. To 
construct a solution in a given domain it is required to 
specify (i) the boundary of the domain; (ii) appropriate 
boundary conditions, e.g. n • B = 0, where n is the unit 
normal; and (iii) homological conditions, i.e. line inte- 
grals (fluxes) around topologically inequivalent loops. To 
specify a solution in a simple torus it is sufficient to spec- 
ify \i and the toroidal flux, while in a doubly-connected 
annulus the poloidal flux must also be specified [43] . 

Solving the Beltrami equation in general C a is, how- 
ever, an intractable numerical problem. Because of the 
topological complexity resulting from the infinity of is- 
lands embedded in the chaotic sea, there is an infinity of 
inequivalent closed loops. The outer boundary of each 
chaotic sea is presumably a critical KAM torus, which is 
not smooth, and the normal to the fractal boundary is 
not defined. 

Furthermore, a continuous, non-trivial pressure, that 
is consistent with a generic non-intcgrablc field, must be 
fractal. To sec this, we may assume that a finite pressure 
gradient is supported by the KAM tori. The Diophantine 



condition serves as a simple, proxy indicator function de- 
scribing the existence of KAM tori in the fractal phase 
space of a generic non-integrablc field (though the more 
complicated Bruno function [44] is probably a better ap- 
proximation). Let us consider a Diophantine pressure 
profile, p(i), defined p'{t) = 1 if \t — n/m\ > r/m k for all 
integers n and m, and p'(t) = otherwise, supplemented 
with the condition p(0) = 0. The function p(i) is con- 
tinuous by construction (i.e. the derivative is nowhere 
infinite) and we assume that r and k have been chosen 
so that p'{i) is non-zero on a set of finite measure (so 
that not all the excluded regions overlap) so that p(i) is 
non-trivial. Even for this 'toy' model, numerically ap- 
proximating the function p(t) given p'(t) is rather com- 
plicated. 

An approximation to p(t) may (in the case of contin- 
uous p 1 ) be constructed using a tagged partition, i.e. 

~ HiP'{ x i)( l i - where x t G [t;-i,ii] and 

= to < *i < • • • < *jv = *• However, because p'{t) 
is nowhere continuous except where p'(t) = 0, the result 
depends on the choice of x% even when |*j — — > 0, Vi. 
That is, the Ricmann integral of p'(t) does not exist. The 
error between this approximation and the exact solution 
is not bounded. 

More sophisticated numerical discretizations could be 
derived; for example, by choosing the Xj in the tagged 
partition to coincide with the locally most irrational, 
i.e. by constructing an 'irrational' tagged partition [45]. 
However, a precise treatment would involve numerically 
approximating the Lebesgue integral and complicated 
measure theory. In 1967, Grad [46] made a similar com- 
ment, describing the pressure as "pathological" . 

Furthermore, the fractal structure of non-integrable 
fields in toroidal confinement devices will be far more 
complicated than that described by the Diophantine con- 
dition, and may not (will not!) be known apriori. There 
are numerical diagnostics for determining the structure of 
phase space, such as Greene's residue criterion, but these 
diagnostics come at considerable computational cost. 

Considering (i) that a nonlinear equilibrium calcu- 
lation will inevitably require an iterative approach, in 
which the fractal phase space structure of the field may 
need to be re-evaluated at each iteration, and (ii) that the 
critical KAM tori are fragile, and an infinitesimal change 
in B can cause an abrupt, finite change in the volume of 
any given chaotic sea, and (iii) that the fractal structure 
of phase space will need to be resolved sufficiently accu- 
rately in order to guarantee that an appropriately defined 
error is below some bound, that can be made arbitrarily 
small as the numerical resolution is increased; we may 
expect that this computational cost would be excessive. 

For our purpose of constructing a robust and efficient 
numerical solution of well-defined, 3D MHD equilibria 
with non-integrable fields, with a bounded error that can 
be made arbitrarily small, it is far better to work with 
smooth functions, and to employ an algorithm that does 
not depend on resolving the infinitely complicated struc- 
ture of phase space. So, we extend our class of functions 
for p and B beyond globally continuous functions, as non- 
trivial, continuous functions that satisfy force balance are 
necessarily fractal, and consider instead functions that 
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are continuous, and smooth, almost everywhere; that is, 
we consider functions that are smooth except for a fi- 
nite set of discontinuities, which can be easily managed 
numerically. 



B. weak solution on discrete partition 

Continuous pressure profiles are not the most gen- 
eral solutions of Eq. (1). Discontinuous pressure pro- 
files may seem unphysical, but they are a valid solution 
class within the zero-Larmor-radius MHD model we have 
adopted. If continuous, globally smooth solutions are re- 
quired, then additional 'non-ideal' physics should be in- 
cluded [47]. For example, including a small, but finite, 
diffusion of the pressure perpendicular to the magnetic 
field will provide solutions with a globally smooth pres- 
sure; and including a small resistivity will prevent the 
formation of singular currents. Appropriate source terms 
are required to balance dissipativc effects. 

This is the approach adopted by various codes [48-51] 
that can approximate an MHD equilibrium as a resistive 
steady state, but which are best described as initial-value, 
time-evolution codes and cannot, strictly, compute an 
equilibrium that satisfies B • Vp = 0, with the pressure 
given. The algorithms these codes employ become in- 
creasingly ill-conditioned as the non-ideal terms approach 
zero [52]. 

We now describe a restriction of the solution class that 
greatly simplifies the equilibrium problem. A discrete 
invariant partition of phase space is constructed. The 
disjoint, invariant sets that are surrounded by a given 
primary chaotic sea, e.g. the hierarchy of island chains, 
are absorbed into the chaotic sea itself, which then be- 
comes cither a simply or doubly connected region; and 
the outer boundary of these regions is extended past 
the adjacent, critical boundary surface to a smooth, no- 
ble surface, which also serves as the boundary for the 
next 'extended' chaotic sea. That is, we choose a set of 
smooth, noble KAM tori, X;, where 2 = 1,2,... Ny, that 
partitions phase space into Ny invariant toroidal or an- 
nular subvolumcs, V/. Each Vi is an invariant set under 
the field line map, but not necessarily an ergodic invari- 
ant set because the field may not be totally chaotic. 

In each region, V;, we equate all the (j, a > to a single 
constant fii, and all the p a ' and ppi to a single constant 
Pi , where a' and /3' label all the chaotic seas and invariant 
tori within V;. Each V/ is simply or doubly connected 
with a smooth boundary and it is a simple computational 
task to solve VxB| = in each V/. We will enforce 
the constraint that n • B = on the 1/ , but otherwise the 
topology of the field in each V/ is unconstrained. 

For the pressure, rather than restricting attention to a 
globally continuous pressure with finite pressure- gradient 
on the uncountably infinite Cp, we instead consider a 
piecewise continuous pressure with finite pressure- jumps 
on the finite set X/. Intuitively, wc imagine that all of 
the pressure in the continuous-but- fractal model that is 
supported by the Cp in the vicinity of a selected noble 
KAM torus is placed on the noble torus itself: all of the 
pressure is placed on a finite selection of the most irra- 



tional surfaces. The vanishing of the divergence of the 
stress tensor, Eq. (1), in a neighborhood of a surface of 
discontinuity gives a condition [1] that must be satisfied 
at the interfaces, namely that the total pressure must be 
continuous across the X/, i.e. \[p + B 2 /2]] = 0. 

We have described the interfaces where there exists 
a discontinuity in the pressure and the tangential field 
as KAM surfaces, but this is rather loose terminology. 
Such interfaces are perhaps 'double-sided' KAM surfaces, 
being covered by an field line = integral curve with ir- 
rational frequency, of the field, B„, immediately inside 
the torus, while also being covered by an integral curve 
with the same irrational frequency of the perhaps dif- 
ferent field, B + , immediately outside the torus. In the 
next section, where we describe the MRXMHD energy 
functional, wc shall refer to the X; as ideal interfaces and 
the Vi as relaxed volumes, and describe why the Xj are 
required to have irrational rotational transform. 

In the above discussion, we have argued that stepped- 
pressure equilibria arise naturally when one seeks a nu- 
merically tractable discretization of the equilibrium prob- 
lem that is consistent with the zero-Larmor-radius model 
of MHD; satisfies B • Vp = 0; and is consistent with what 
is known about the fractal phase space structure of non- 
integrable fields. The equilibria could also be described 
as multi-volume sharp-boundary states. A major motiva- 
tion for pursuing this model is that Bruno & Laurence [4] 
have proven that such stepped-pressure equilibria exist 
(provided the departure from axisymmetry is sufficiently 
small). The number of volumes, Ny, and interfaces may 
be made arbitrarily large. We can, depending on the nu- 
merical resources available, consider a sequence of invari- 
ant partitions with increasing Ny, so the discontinuities 
in the pressure are made arbitrarily small, in order to 
study the nature of a continuous-but-fractal equilibrium 
via a sequence of well-defined, stepped-pressure equilib- 
ria. Stepped-pressure profiles are sufficiently general to 
represent observed profiles to within experimental error. 

In order to explore the properties of these equilibria 
in arbitrary geometry, i.e. to go beyond what may be 
proved analytically, this model has been implemented 
numerically in the stepped-prcssure equilibrium code, 
SPEC, as will be described below in Sec. IV. We now 
show that there is an multi-region, relaxed MHD energy 
functional, that we call MRXMHD, that has stepped- 
pressure equilibria as extremizing solutions. 

III. ENERGY FUNCTIONAL METHOD 

The classic MHD energy functional [13] is given by the 
integral 

where V is the plasma volume bounded by a toroidal 
surface, dV. Ideal equilibria arc obtained when the 
plasma is in a minimum energy state: more pre- 
cisely, when the energy functional is extremized allow- 
ing for a restricted class of variations, namely ideal 
variations. The equation of state, dtip/p 1 ) = 0, where 
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d t = dt + v ■ V and v is the 'velocity' of an assumed 
plasma displacement, v = d t £, may be combined with 
mass conservation, dtp + V • (p v) = 0, to obtain an 
equation that constrains the variation in the pres- 
sure, Sp = (7 — 1)£ • Vp — 7V • {p£). Faraday's law, 
<9tB = V x E, may be combined with the ideal Ohm's 
law, E + v x B = 0, where E is the electric field, to ob- 
tain an equation that constrains the variation in the mag- 
netic field, <5B = V x (£ x B). Note that this last con- 
straint does not allow the topology of the field to change. 

The first variation in the energy due to an ideal dis- 
placement, £, that is assumed to vanish on the boundary, 
is given by 

SW= I (Vp-jx B)-£dv. (7) 
Jv 

Extremizing solutions satisfy the ideal force-balance con- 
dition, Vp = j x B. In order to uniquely define an equi- 
librium, in addition to the shape of the plasma bound- 
ary, it is required to specify the pressure, and either 
the rotational-transform or the parallel current density 
[9, 53]. (Note that the constraints of ideal MHD places a 
constraint on the differential toroidal and poloidal fluxes 
and the rotational-transform, namely dipp/dipt = *•) 

The VMEC [54] and BETAS/NSTAB codes [55, 56] 
are based on this approach. These codes assume that the 
magnetic field is intcgrablc and allow for smooth pressure 
and rotational-transform profiles. 

In general 3D geometry, there is a singularity in the 
resonant harmonic of the parallel current in equilibria 
with nested flux surfaces [57]. Writing the current as 
j = ctB + the quasineutrality condition, V • j = 0, re- 
quires that the parallel current must satisfy the mag- 
netic differential equation, B • Vtr = —V • jx, where we 
may consider the perpendicular current to be driven by 
the pressure gradient, ji = Bx \7p/B 2 . Magnetic dif- 
ferential equations are densely singular [58]. The sin- 
gularity may be exposed, in the integrable case, by 
the use of straight field line coordinates, which allow 
the directional derivative along the magnetic field to 
be written y^B • V = idg + d^. Using a Fourier repre- 
sentation, e.g. a = J2 m n a m,n exp(im8 — in<f>), we de- 
rive a m>n = -i(y/g V • j±) m , n /(mt - n)+j m , n 8(rm- n). 
The first term is called the Pfirsch-Schliiter current and 
has al/i style singularity at the rational surface, where 
x = t — n/m. The second term, the 5-function current, 
is generally required to 'shield' out resonant magnetic 
fields that would otherwise destroy the nested family of 
flux surfaces (a more precise discussion of the ^-function 
current is provided in [59]). 

In general geometry, the only way to avoid the 1/x 
singular currents is ensure that the pressure gradient is 
zero in the vicinity of the rational surfaces (or to ensure 
that no rational surfaces are present). As the rational 
surfaces arc dense in space, to avoid the 1/x singularities 
the pressure gradient must be zero everywhere. 

(Despite these concerns near the rational surfaces, 
VMEC, in particular, does an impressive job of robustly 
constructing global approximations to 3D equilibria with 
arbitrary pressure profiles; presumably, this is because 
VMEC seeks approximations to minima of the global en- 



ergy functional, and does not directly seek solutions to 
= j x B pointwise.) 

A. MRXMHD energy principle 

The first step towards constructing the multi-region, 
relaxed MHD energy functional is to partition space 
into discrete volumes. We introduce a set of nested, 
toroidal surfaces, Ii, for 1 = 1,2,... Ny where Xi = dV 
for I = Ny. The energy local to each volume is 

where Vi is the toroid enclosed by Ii , and V; is the an- 
nular volume enclosed by and X\ for 1 = 2,... Ny. 

We again assume that the plasma is in a minimum en- 
ergy state; however, we allow for the effects of small resis- 
tivity: in each Vi , the magnetic field may relax and recon- 
nect (and so topological constraints between the toroidal 
and poloidal fluxes, the rotational-transform, and the he- 
licity are broken). But, in order to retain some control 
over the equilibria, we consider the I; to be preserved as 
ideal barriers that restrict both pressure transport and 
field transport. Rather than continuously constraining 
the topology, the topology is discretely constrained. This, 
or something equivalent, is required in order to avoid triv- 
ial solutions. 

In each Vi, the mass and entropy constraints usually 
used in ideal MHD do not apply to individual fluid ele- 
ments, but apply instead to the entire volume, giving the 
isentropic, ideal-gas constraint, 

PiVf = a u (9) 

where Vi is the volume of V; and ai is a con- 
stant. The internal energy in V; is J v pi/{"i — 1) dv = 

aiV, _7 V(7 ~~ l)i ail d the first variation of this due to a 
deformation, of the boundary is —p J gv £ ■ ds. 

To constrain the relaxation of the magnetic field in 
each V/, we follow Taylor [60], who argued that the 'most 
conserved' invariant for a weakly resistive plasma is the 
helicity [5, 61], 

Ki = [ A B dv, (10) 
Jv, 

where A is a vector potential, B = V x A, which we 
consider to be diffcrentiable and a single- valued function 
of position. The helicity is related to the Gauss linking 
number: it reflects how 'knotted', or 'twisted', the mag- 
netic field lines are [5, 42, 62]. The helicity in Eq. (10) 
is not gauge-invariant. A gauge-invariant form is con- 
structed by adding the loop integrals Aip p <f> s A ■ dl and 
Aipt §r -A- ' dl, where S is a poloidal loop on and C 
is a toroidal loop on I; . 

In each V;, variations in the pressure and the field, 
and the geometry of the interfaces, are allowed in or- 
der to extremize the energy functional. These variations 
are arbitrary, except for (i) the mass-entropy constraint, 
PiVp = const; (ii) helicity conservation in each V/; (iii) 
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the interfaces must remain tangential to the magnetic 
field; and (iv) the magnetic fluxes are conserved. 
The MRXMHD energy functional is 

F = J2[Wi-^(K l -K\ )]. (11) 
I 

The hclicity constraint, Ki = Ki_ where Ki i0 is a given 
constant, is enforced explicitly by introducing a Lagrange 
multiplier, in each V/. The flux constraints and the 
tangentiality condition at the interfaces will be enforced 
implicitly by constraining the representation of the mag- 
netic field. 

The most general function space for B in each volume 
is space of vector-valued functions whose magnitude is 
square integrable, i.e. B £ C 2 (Vi), by which it is meant 
that B 2 £ £ (V;), and jC 1 (Vi ) is the standard notation 
for the space of integrable scalar functions. More pre- 
cisely, we follow Yoshida et al. [6, 40] and restrict B to 
C 2 (Vi), which they define as the subspacc of £ 2 (V/) oc- 
cupied by divergence-free fields that obey B • n = on 
the boundary. Similarly, the pressure is required to be 
integrable, i.e. p £ £ 1 (V;). While these are the least re- 
strictive spaces required for a weak formulation, in order 
for the solutions to obey tractable local differential equa- 
tions almost everywhere we will, after deriving the Euler 
Lagrange equations for states that extremize Eq. (11), 
further restrict the allowed function spaces by assuming 
that p is piecewise constant and that B piecewise satis- 
fies a simple elliptic partial differential equation, which 
is solved numerically using a mixed Fourier and finite 
element method. 

The variation in the 'local' constrained energy func- 
tional, Fi = Wi — hi (Ki — Ki t0 ) /2, due to arbitrary vari- 
ations in the field, SH = V x 8 A, and arbitrary varia- 
tions, in the interface geometry, is given by 

SFi = { (VxB- W B)-dAcfo 

JVi 

- f ( Pl + B 2 /2)Z-ds. (12) 

JdVi 

The variation in the magnetic potential, <5A, is free 
within V/ , and so within each V/ the topology of the field 
is arbitrary; but at the Ii it must obey 

nx6A = -n-£B + nx WSg, (13) 

so that n-B = remains satisfied; and where Sg(r) is the 
variation in a single- valued gauge potential, g, required 
for generality but physical quantities are invariant with 
respect to gauge choice. The line integrals of A along 
arbitrary loops £ po1 and £ tor are related to the poloidal 
and toroidal magnetic fluxes. 

The enclosed toroidal fluxes in each volume and the 
poloidal fluxes in each annular region constrain the mag- 
netic field from being trivial. We use gauge freedom to 
specify the loop integrals of § A-dl on each interface. For 
gauges satisfying these conditions, the difference between 
the gauge-invariant hclicity and the gauge-dependent he- 
licity is a constant. 

The Euler-Lagrange equation for F to be stationary 
with respect to variations in the magnetic field in each 



V/ is the Beltrami equation, V x B = /J;B. The Euler- 
Lagrange equation for F to be stationary with respect 
to variations in the interface geometry is that the to- 
tal pressure must be continuous across the interfaces, 
\[p + B 2 /2}} = 0. States that extremize the MRXMHD 
energy functional arc stepped-prcssure equilibria. 

The pressure and tangential field are discontinuous 
across the interfaces, but these comprise a finite set of 
measure zero and so p and B 2 are both integrable func- 
tions: the model is consistent with our goal of construct- 
ing weak solutions via an energy-integral approach. The 
discontinuities are easily accommodated for in the nu- 
merical discretization; within each volume a continuous, 
smooth representation for the vector potential is allowed. 

To avoid a problem with 'small denominators', as will 
be discussed below, we will typically enforce the condi- 
tion that the interfaces have irrational transform. The 
problematic Pfirsch-Schlutcr currents are eliminated be- 
cause the pressure gradient is identically zero across the 
resonances. The <5-function currents are also not present 
because the topology of field is unrestricted at the ratio- 
nal surfaces, i.e. magnetic islands are allowed to form. 

There are, instead, a finite set of surface currents at 
the irrational, ideal interfaces, given by k = [[B]] x n, 
where [[B]] is the tangential discontinuity in the field. 
These are required to enforce the topological constraint, 
the topological constraint in this case being that a noble 
irrational surface exists; and the topological constraint 
is required so that the magnetic field matches the given, 
stepped-pressure profile. Given the KAM theorem, this 
topological constraint is presumably easier to enforce [63] 
than forcing a rational flux surface. We expect these cur- 
rents to be dominated by the discontinuity in p, and may 
be thought of as a discrete approximation to the pressure 
induced currents. These irrational surface currents may 
be compared to (but are different from) the (5-function 
currents shielding at the rational surfaces, which are re- 
quired in the linearly-perturbed, ideal-equilibrium codes 
IPEC [64] and CAS3D [65]. The (5-function currents at 
the rational surface currents do not describe the pressure 
driven l/x singularities. 

We invoke multi-region energy minimization with he- 
licity conservation primarily mathematical device 

to achieve a variational formulation of the restricted 
equilibrium class, namely stepped-pressure equilibria. 
MRXMHD is, however, a generalization of the variational 
principle enunciated by Woltjer [66] to generate linear 
force-free fields of interest in astrophysics, and developed 
by Taylor [60] to model fusion plasma experiments. 

The success of the Taylor relaxation theory in describ- 
ing experimental data suggests that the MRXMHD ap- 
proach may likewise aid physical interpretation of partial 
relaxation, reconncction and self-organization in toroidal 
plasmas supporting a non-trivial pressure profile. Unlike 
Taylor's globally relaxed model, which gives a constant 
pressure across the plasma, MRXMHD is only locally re- 
laxed, i.e. it is partially constrained; arbitrarily many 
interfaces may be included, each with an associated ideal 
= topological constraint. 
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B. transform constraint 

A close examination of the interface force-balance 
condition, [[p + B 2 /2]} = 0, reveals a Hamiltonian 
system, which we call the pressure-jump Hamilto- 
nian [1]. Let p_ and B_ be the pressure and field 
immediately inside a given interface and p + and B + 
be the pressure and field immediately outside. By 
combining (i) the general, covariant representation 
for the field, B = B s \7s + B g V9 + B^Vcj), with (ii) 
V x B = /iB, and (iii) the tangentiality condition, 
B • n = 0; we may write Be = defe and B^, = d^f , and 
B 2 ={g<j } <j } fefe - ^ge^fefcf, + geef<pf<j>)/(geeg<j><j> - g6<i,g8<p), 
where /(#,</>) is a surface potential and gee, ge^ and 
g^cf, are metric elements (local to the interface). Now, 
consider the case where both B and the geometry of 
the interface are known, and we seek a solution for B + 
that satisfies H = const, where 



H = 2{p_-p + ) = B\-B 2 _. 



(14) 



We may write H = K{9, </>, fe, f<p) + V(6, <j)), where 
V = —B 2 _ is assumed known, and fe = def and 
f<t> = d^f, where / is an as-yet-unknown surface potential 
for B+. 

To derive the solvability condition we treat fe and 
as independent quantities (generalized momenta) 
and recognize H as a 2-dof Hamiltonian with a 
conserved energy, 2(p_ — p+). Then, H = const 
along a trajectory given by Hamilton's equations: 
de/dt = dH/dfg, dfg/dt = -dH/d6, defy/dt = OH/df^, 
and df^/dt = —dH/d(f>; where t is an artificial 'time'. 
This system may be reduced to a 1^-dof system by using 
4> as the time-like integration parameter (always possi- 
ble if dt4> 7^ 0) and eliminating the integration of f$ in 
favor of inverting K(9, (f>, fe,f<t>) + V(6, <p) = 2(p_ - p + ) 
for f^, so that is assumed to be a function of 9, <fi 
and fe, i.e. = f$(9, <fi, fe), where the dependence on 
2(p_ — p+) is implicit. The trajectory is then described 
by 9 = dt9/d t 4> and fe = dtfe/dt4>, which may be inte- 
grated in <f) from an initial starting point, (9,fe), on 
a Poincare section, e.g. <fi = 0. If the trajectory lies 
on an invariant surface then it is possible to construct 
fe = /e (#,</>) j and fe and recover their interpretation 
as derivatives of a surface function: there exists a well 
defined f(9, </>) such that fe = def and = d^f . That 
is, if the trajectory lies on an invariant surface, then a so- 
lution for B + that satisfies 2(p_ — p + ) = B\ — B 2 _ may 
be constructed. 

An invariant surface can only exist if it avoids the prob- 
lem of small divisors. Note that 9 = B 9 / B^ , so there is 
a fundamental relationship between the pressure-jump 
Hamiltonian and the field-line Hamiltonian; and force- 
balance can only be satisfied if the rotational transform 
of the interfaces is irrational [1]. 

Many authors have considered sharp boundary equilib- 
ria either theoretically or in simplified geometry [2, 43, 
67-77]. This paper, and our earlier paper [78], represents 
the first numerical study of toroidal 3D equilibria with 
multiple Beltrami regions within the plasma. The only 
3D calculation of which we are aware is the early paper 
of Betancourt and Garabedian [79] , who consider a free- 



boundary problem with both the vacuum region and the 
plasma being Beltrami regions with /i = 0. 

A number of 3D MHD equilibrium codes based on the 
assumption of continuity and differentiability of p and B 
have been written [47, 53, 54, 56, 79-89]. These have 
cither constrained the magnetic field to be globally in- 
tcgrable; have not employed numerical algorithms that 
explicitly accommodate the singularities in the parallel 
current at the rational surfaces; do not constrain the 
profiles, and allow the pressure, current and transform 
profiles to 'evolve' during the calculation in a fashion 
more akin to initial-value, time-evolution codes rather 
than what is suitable for an equilibrium code; have in- 
troduced small non-ideal terms, so that the B • Vp 7^ 0; 
have ignored the fractal hierarchy of the crgodic invariant 
sets; or employ ill-posed numerical algorithms (e.g. the 
so-called Spitzcr [18] iterative approach [85, 90], which 
attempts to invert densely-singular magnetic differential 
equations [91]). While they have produced a variety of 
results, their lack of formal foundations leads them to 
fall short of the numerical rigor (e.g. demonstration of 
convergence, quantification of error, estimate of stability 
[92]) available in the axisymmetric case. 

We will now describe the stepped-pressure equilibrium 
code, SPEC, and demonstrate that the solutions are well 
defined by presenting convergence studies. 



IV. NUMERICAL DISCRETIZATION 

A Fourier representation is employed for all 
doubly-periodic, scalar functions. Even functions, 
f{-9,-Q) = f{9,Q), are written 

N 



f = ^2fo.nCOs(-nN P () 

M N 

+ Y. Y. fm,nCos{m9-nN P 0, (15) 

m— 1 n—~N 

where Np is the field periodicity. The resolution 
of the Fourier representation is determined by M 
and N, and the total number of Fourier harmon- 
ics is N M n = (N + 1) + M(2N+ 1). The poloidal, 
9, and toroidal, £, angles are, as yet, arbitrary. 
The Fourier summation will be written concisely as 
/ = J2j fj cos(nij9 — nj(), where (mi, n{) = (0, 0), etc. 
A similar description is used for odd (i.e. sine) functions, 

f(-o,-o = -mo- 

An initial guess for the geometry of a set of Ny nested, 
toroidal surfaces, I/, is assumed given. For expedience, 
we assume stellarator symmetry [93] , so that the Xi may 
be described by R(9, QR + Z(9, () k, with 



3 



(16) 



where R = cos <f> i + sin 4> j, and i, j and k are the Carte- 
sian unit vectors. 
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To enforce various boundary conditions, it is con- 
venient to use toroidal coordinates, (s, 0, (,), that are 
adapted to the interfaces. These coordinates are defined 
inversely via 

R = R(s, 6,(), <t> = -Q, Z = Z(s, 0, C). (17) 

The Jacobian of the (s,9,() coordinates is 
^fg = R(R s Ze — RgZ s ). The 'lower' metric coeffi- 
cients, g a fi, arc given by g a p = R a Rp + Z a Z/3 + 5 a pR 2 , 
where S a p = lifa = /3 = £ and 8 a p = otherwise. The 
coordinate functions are given by 

R(s,9X) = ^ Rjjs) cos(nij9 — iijC,), 

4^ . ( 18 ) 

Z(s,9,0 = y Zj[s)sm.(mj6 — rijQ, 

3 

where the Rj(s), Zj(s) are a piecewisc- linear interpola- 
tion of the i?zj and Zij. (A piecewise-cubic interpolation 
would give a continuous Jacobian across the interfaces, 
but this is not required.) If the toroidal flux is mono- 
tonic increasing, then s = ipt> normalized to its value at 
the outermost interface, is a suitable radial coordinate. 
In this case, s ~ r 2 , where r is a polar-like radial coor- 
dinate. More generally, we may use the interface label 
itself as the radial coordinate, i.e. si = l/Ny- 

In the innermost volume, regularization factors 
must be included to prevent the interpolated co- 
ordinate surfaces from overlapping. These factors 
may be derived by considering an arbitrary, regular 
(infinitely-differcntiable) function, h(x,y), at the origin, 
h(x, y) = h + xh x + yh y + \ {x 2 h xx + 2xyh xy + y 2 h yy ) + 
.... By constructing a Fourier representation, 
H r > °) = Em iKi( r ) cos(m6) + h s m (r) sin(m0)] where 
x = rcos9 and y~rsm9, we obtain after repeated 
applications of double-angle formulae, 

h m (r) = r m {a r° + a x r 2 + a 2 r 4 + a 3 r 6 + ...). (19) 

So, in the innermost volume, s < si, we write 
Rj(s) = Rj,is mj / 2 1 s" 1 ^ 2 , and similarly for Zj(s), where 
s ~ r . 

In Vi that is bounded by and Ii, a general co- 
variant representation of the magnetic vector potential 
is 

A l = A s , l Vs + Ae,iV9 + A u VC. (20) 
To this add a gauge term, Vgi(s, 9, £), where gi satisfies 

d sgi {s,e,Q = -1^(3,6,0, 

deafa-iAC) = -4^(^-1,0,0 + ^-1, (21) 

<9 C g;(s;_l,0,C) = -A fii (si_i,0,C) +1pp,l-l, 

for arbitrary constants, ipt,l-i an d Vp,i— l- Then 
A; = A; + Vgi is given by A/ = Ag,jV9 + A^jVC with 

A e ,i(si-i,0,C) = ipt,i-i, , 00 x 
A u ( Sl ^,Q,Q = ^ppj—i . 

For stcllarator symmetric equilibria, Agj and Aq^ may 
be represented by cosine series, 

A e j(s,9,() = Ag t i ij (s) cos(mj9 - n 3 <), 

j (23) 
Ac,i(s,6,() = ^2 A (,ij( s ) C0S ( m 3 e - n j0> 
j 



where As : ij(s) and A^ij(s) arc represented using finite- 
elements, as will be described below. 
The toroidal flux is given by 

/fids = I A ■ dl = 27r ipt,l-u (24) 

JS JdS 

where the surface S is that part of the ( = plane 
bounded by s = s;_i. The poloidal flux is given by 

/ B ds = <£ A-dl = 27r W/-i- (25) 

JS JdS 

where the surface 5 is that part of the = plane 
bounded by s = s;_i. 

The boundary condition that the inner interface is a 
flux surface becomes B • Vs = 0, which implies 

- m j A C,u( s i-i) - njAg M {ai-x) = 0. (26) 
Combing Eq. (22) and Eq. (26) we have 

^w(-w) = {V ; )ll (27) 

A ( ,A^i)={\ 1 ; (28) 

The condition that the outer interface is a flux surface is 
satisfied by writing 

A e ,i(s l )=d e fi(6,0, A a ( Sl ) = d c M9,0, (29) 
for arbitrary / of the form 

fl = 4't,i8 + ipp.iC + ^2 fl,j sm(m 3 -0 - TijC), (30) 
o 

and ip t< l = Ag,i,i[si) and ip p j = Aa (j i(s;). We have 

A >" M ={ii : j -I: (31 » 
A «^={-t% : i ■!: (32) 

The radial dependence of the vector potential 
harmonics is described using finite-elements. A 
continuous function, f(x), with x€ [0,1], may 
be approximated using the linear basis functions, 
1 Pl,o{x) = 1 — x and ipnfi(x) = x, according to 

f(x) = fL,0<fL,o( x ) + fR,0<PR,o( X )> where fL.O = /(0) 

and fufi = /(l). A piecewise-linear interpolation of the 
vector potential gives a discontinuous magnetic field in 
each V/. While this is legitimate as far as the energy 
integral is concerned (the magnetic field remains an 
integrable function), we prefer a smoother interpolation. 

For piecewise-cubic interpolation, the basis functions 
are (fL,o(x) = 2x 3 — 3x 2 + 1 and (fL,i(x) = x 3 — 2x 2 + x, 
and their 'reflections', ipn tP (x) = (— 1) p <^l )P (1 — x). An 
arbitrary smooth continuous function is approximated 

by fi x ) = J2p"o\-fL,P ( PL,p{ x ) + fR,p<PR,p(®)], where 
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f L:1 = f(0) and f R>1 = f'(l) and N D = 1. For 
piecewise-quintic interpolation the same expression ap- 
plies, but with Ne> = 2 and >pl.o( x ) = — 6x 5 + 15x 4 — 
10x 3 + l, = -3x 5 + 8x 4 -6x 3 + l, and <y9i, 2 (s) = 

— \x 5 + |x 4 — §x 3 + \x 2 , and their reflections. 
In each Vz, a regular, radial sub-grid is established, 

si,i = 8i-! +i(si -si-x)/Ni, (33) 

for i = 0, 1, . . . , iVj. The resolution, A;, of the radial sub- 
grid may be different in each Vi- The vector potential 
harmonics are written 

Agjj(s) = 2J [^Vzj, P ,i-i <Pl, p (x) + ¥r, p { x )\ > 

where x = (s - s; ; i_i)/Asi with As; = sij - sj^-i. 

The vector potential is completely specified by Ag t ij, Pt i 
and ^4^./j.p.i, which are the p-th. derivatives of the 
(m,j,rij) harmonics of Ag and Aq on the i-th grid-point 
in the i-th annulus. Except for the subtlety required to 
ensure the field is tangential to the outer interface, see 
Eq. (31) and Eq. (32), these are the independent pa- 
rameters that describe the vector potential, and thus the 
magnetic field. 

In the innermost volume, the condition that the field 
is tangential to the inner interface is replaced (because 
there is no inner interface) by the condition that the vec- 
tor potential is analytic at the coordinate origin. Assum- 
ing s r~j r 2 , we may enforce regularity at the origin and 
restrict the gauge by including s m i/ 2 radial factors with 
the Agj j tP i and ij Pt i, with the boundary conditions 

Ae,!,j,o,o = for all j, (34) 
A^x,j,o,o = for uij = and rij 7^ 0. (35) 

The mixed finite-element, Fourier representation of the 
magnetic vector potential is inserted into the MRXMHD 
energy integral, F = JZj-Fj, where the 'local' energy 
functional is given by Fi = Wi — m (Ki — Ki Q ) /2, where 
Wi = J [p/(7- 1) + B 2 /2_ \dv and K t = jA-Bdv. 
With A = Ag\79 + A^\/(, the magnetic field 
B = B s e s + B e e g + B^e c is 

v^B = {d 6 A c - d c Ag) e s - d s A ( e e + d s A e e c , (36) 

and 

B 2 = B S B S g ss +2B s B e g s g + 2B s B^ g s( 
+ B e B e g ge + 2B 9 B< g ec + B<B< g co 
V?A-B = -Agd s A c + A c d s Ag. 

After substituting these into the local energy and helic- 
ity integrals, the Ag.ij. p .i and A^j^^i are multiplied by 
terms such as 

f'% fd9 [£<plA*) <PR, q (*)fMC), (37) 

Jsi,i-i JO JO 

where, for example, / = sm(rrij9 — rijQ) g a p cos(mkO — 
UkC)t an d g a p are the metric elements (note that the 



g a p depend on the Rij and Zij that define the inter- 
face geometry). These integrals arc computed by con- 
structing a fast Fourier transform of f(s, 9, £). The inte- 
grals over the angles then become trivial and we obtain 
f(s) = (pL, P (s)(pR,q{s) Jd9j d(f(s,6,(). The remaining 
radial quadrature is approximated using Gaussian inte- 
gration, 

/ (38) 

where the 'weights', Wi and the 'abscissae', S4, are cho- 
sen to optimize accuracy, and Nq is a numerical resolu- 
tion parameter that depends on the order of the poly- 
nomial being integrated, i.e. Nq depends on the order 
of the finite-element basis expansion for the Agjj tPj i and 
^C,i,j,p>£) and the order of the coordinate interpolation of 
the Rij and Zij. 

In each V;, the local energy functional, 
Fi = Wi — m (Ki — Ki t0 ) /2, depends on the interface 
geometry, the vector potential, and various input param- 
eters. Specifically, each Fi depends on x = {Ri,j,Zi,j}; 
the enclosed toroidal flux, Aip t .i = V 1 *,; — ipt,l-i] the en- 
closed poloidal flux, Aip p j = ip p j — tpp,i-i (except in the 
innermost volume); the required helicity, -fTz, ; the vector 
potential, aj = {Ag t i^ tPji , A^j jPti }; and the Lagrange 
multiplier, /i;. We may indicate the dependence of F as 

Fi =i^AV> M) A^ Pii ,i^ ,x; W ,a(]. (39) 

In the MRXMHD model, the enclosed fluxes and the he- 
licity in each V/ are assumed given, i.e. the Aip t ,i, Aip p j 
and Ki o are required input parameters. (Note: if the 
interface rotational-transform constraint is to be given 
priority over the conservation of poloidal flux and helic- 
ity, then A-0 Pj ; and A; must generally be allowed to 
vary.) The computational task is to then find extrema 
of F = J2i Pi with respect to the interface geometry, 
x = {Ri,j, and the Lagrange multipliers, /!;, and 

the vector potentials, {Agj j tPl i, A^i ,j, p ,i}- 

The basic algorithm is to consider /i; and 
a; = {Ag t i t j tPt i, j4f,i,j,p,i} to be functions of the in- 
terface geometry and the A%p t ,i, Aip Pi i, and Ki M . That 
is, first, the Beltrami field, V x B = /i/B, in each Vj is 
constructed. We then may write Ft = Fi[Aip p j, A'/ jD ,x], 
where the dependence on Aip t ,i is implicit. (Later, for 
computational expedience, we shall modify this slightly 
by using //; to parametrize the solutions to the Beltrami 
fields and remove AT; j0 , so that F; = Fi[Ax/j Pt i, /x;,x].) 
Then, 'global' equilibrium states are then constructed 
by extremizing F = Fi with respect to the interface 
geometry, x = {R Lj , Z Lj }. 



A. solving VxB = fiB for B, given geometry 

Assuming that the geometry of the interfaces, x, is 
given, there are various numerical methods that may be 
employed to construct the extremizing fields. The first 
method is the standard Lagrange multiplier approach: 
a multi-dimensional Newton method is used to find an 
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extremum of the local constrained energy functional. The 
solution satisfies 



m 

dai 



= 



dm 



(40) 



where, in addition to a; = {-Ae,i,j,p,i, ^-C.,l,j,p,i\i t ne La- 
grange multiplier is explicitly treated as an independent 
degree of freedom and must be adjusted to satisfy the 
helicity constraint. 

This approach cannot distinguish states that minimize 
Wi from states which are saddle points or local maxima 
of Wi. When bifurcated solutions exist, i.e. when there 
exist multiple stationary points of Wj and hence Fi, a 
gradient-descent algorithm such as sequential quadratic 
programming [94] may instead be used to ensure that 
the constructed solution is strictly a local minimum of 
Wi subject to the constraint Ki = K^ Q 

Another method for constructing the Beltrami fields, 
on which we hereafter concentrate, is to assume each B; is 
parametrized by the enclosed fluxes and \ii. The required 
value for the helicity, K^ , may be dropped from the local 
energy functional, to obtain Fi = Wi — mKi/2, and we 
write 



Ft = Fi[Aijj Pt i,fj,i,x.;aLi] 



(41) 



where the dependence on Aip t ,i is implicit. 

The local energy functional, Fi, is quadratic in the 
Ae,ij, P ,i and A{j t j >Pt i, and the 'local' equilibrium condi- 
tion, dFi/d&i = 0, gives a system of linear equations to 
be solved for the vector potential. We call this the Bel- 
trami linear system, as it is analogous to V x B = /i;B 
and can be represented as 



G • a = c, 



(42) 



where the matrix G depends on the geometry, the fluxes 
and Hi, i.e. G = G [A^ p ,z , A*z , x] ; and similarly for the 
right-hand-side vector, c. This system is inhomogeneous 
(i.e. c is non-trivial) because of the 'driving' terms 
-4fl,i,/,o,o = i>t,l and A^i.i.o.o = VW- Given x, there is 
a two-dimensional family of solutions; each solution is 
parametrized by Aip p ^ and /i;. 

There is an abundance of numerical methods and 
'canned' numerical routines available for solving linear 
equations, and any mathematical structure present can 
be exploited. For example, usually the matrix G is posi- 
tive definite and it is typically very sparse, and an initial 
'guess' for the solution is often available (particularly so 
during an iterative calculation). Employing numerical 
methods that exploit the sparsity and positive definitc- 
ness can significantly improve code performance. 

It will be efficient to know how the Beltrami field in a 
given volume varies with small variations in both the in- 
put parameters and the interface geometry. In equilibria 
that are globally constrained by ideal MHD, variations 
in the magnetic field are related to variations in geom- 
etry via <5B = V x (£ x B). In our case, we can com- 
pute the change in the vector potential resulting from 
an infinitesimal change in Aip p j, or in /!;, or in the in- 
terface geometry, x, using matrix perturbation theory, 
(G + dG) ■ (a + da) = (c + dc), so that to lowest order 



G • da = dc — dG ■ a. 



(43) 



The infinitesimal variations, dG and dc, resulting from 
infinitesimal variations in Aip^ p and hi are rather simple 
to write down because Aipi tP and \ii just appear as factors 
multiplying various geometric quantities in G and c. The 
derivatives with respect to the Rij and Zij are more 
complicated as these involve differentiating Eq. (37). 

The Beltrami field in V; depends on the geometry of 
both the 'inner' interface, i.e. the Ri-ij and Z\—\ j, and 
the 'outer' interface, i.e. the Rij and Zij; giving a total 
of 4Nmn — 2 geometrical degrees of freedom, where Nmn 
describes the Fourier resolution. In preference over re- 
peatedly inverting the Beltrami matrix ANmn — 2 times, 
which would be the case if finite-differences for example 
were used to compute the change in the Beltrami fields, 
we instead first compute a Cholesky factorization of G, 
i.e. G = LL T . Then, the solution to Eq. (43) is effi- 
ciently given by L • y = b, where b = dc — dG ■ a, and 
L T ■ da = y. 

The helicity, Ki = J A ■ B dv, depends on the solution 
to Eq. (42), which in turn depends on /!;. The helic- 
ity constraint, Ki = K^ , can be enforced by suitably 
adjusting (This is not always possible; this is only 
for configurations in which /i; parametrizes states with 
different helicity.) 



B. transform constraint, noble irrationals 

The rotational-transform constraint can similarly be 
enforced. If only the Beltrami field within a single annu- 
lus is to be constructed, then there is no constraint on the 
allowed values of the transform on the interfaces. Recall 
however, that if the Beltrami fields in multiple volumes 
are to be consistently nested together in a fashion that 
satisfies force balance, an analysis of the pressure-jump 
Hamiltonian derived from [[p+S 2 /2]] = shows that the 
interfaces should have irrational transform. 

We restrict attention to noble irrationals [26], which 
play an important role [12] in the theory of chaos as 
invariant KAM surfaces with noble transform are most 
likely to survive chaos- inducing perturbations [20]. A 
noble irrational is obtained as the limit of an infinite, al- 
ternating path down a Farey tree, which is constructed as 
follows. Begin with a pair of rationals, pi/qi and P2/92, 
which should be neighboring, i.e. \p1q2 — P2Qi \ = 1> and 
without loss of generality we assume that pi/qi < P2/<?2- 
A Farey tree is constructed by successively construct- 
ing the mediants, defined as p/q = (pi +P2)/(qi + 52)- 
This is guaranteed to lie between the original 'parent' 
rationals and so splits the original interval into left, 
[pi/qi,p/q], and right, [p/q,P2/q2], sub- intervals. The 
mediant is neighboring to both parents, and the con- 
struction of the Farey tree proceeds iteratively. An infi- 
nite, alternating path down the Farey tree is a sequence 
of mediants for alternately the left and right subinter- 
vals. Sequences of this type converge to noble irrationals, 
which have continued fraction representations that ter- 
minate in an infinite sequence of l's [12]. It is easy 
to see that alternating paths give Fibonacci sequences 
for the numerator and denominator of the successive ra- 
tionals. For example, beginning from p\/q\ = 0/1 and 
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Vil<li = 1/1 and constructing an alternating path of me- 
diants we obtain the sequence j, j, |, |, |, |, ... 
This allows noble irrationals to be written in the concise 
form b(p ll q 1 ,p 2 ,q2) = (pi +7P2)/(<7i + 792), where the 
golden- mean, 7 = (1 + \/5)/2, is the limiting ratio of suc- 
cessive terms, 7 = F n+ i/F n as n — > 00, of the Fibonacci 
sequence beginning from Fq = and Fq = 1. 

The poloidal angle paramctrization of the interfaces 
is, at this stage, arbitrary; it is not required, and will 
not be required, that the field lines are 'straight'. We 
may, however, construct the straight-ficld-linc angle on 
the interfaces, given the field, by calculating the angle 
transformation, 9 S = 6 + X(9, (), such that 



TABLE I: Flux and transform constraints 



B • V(9 S 
B- VC 



(44) 



where 4 is, as yet, an unknown constant to be deter- 
mined. We restrict attention to angle transformations of 
the form A = -\? sm(rrij9 — nj£), which preserves stel- 
lar ator symmetry but is otherwise general. The Fourier 
resolution of the angle transformation is independent of 
the Fourier resolution, M and N, used to represent the 
interfaces and Beltrami fields, and typically we use an en- 
hanced Fourier resolution for A. With Eq. (36) and using 
B s = 0, Eq. (44) becomes A' e d c \ - A' ( d e \ - A'gt = A' c , 
where the prime indicates radial derivative. By equating 
coefficients, we obtain a system of linear equations for 
the unknowns, A = (*, A2, A3, . . . ) T , which is represented 
as a matrix equation, 



A A = d, 



(45) 



where A and d depend on the A' e and A'^ harmonics 
at the interfaces. Solving this linear-system determines 
A, which gives the rotational-transform on the interface, 
namely t. 

Considering the geometry of the interfaces and the en- 
closed toroidal flux in each V; to be fixed, each Beltrami 
field depends only on /ii and Aip p j. We thus have two 
degrees-of-freedom, and we must satisfy two constraints; 
these constraints being that the field in Vi provides the re- 
quired rotational transform on the inner interface, 
and on the outer interface, X\. (In the innermost vol- 
ume, Vi, there is only one degree of freedom, namely 
jUi, and only one constraint, namely that the field pro- 
vides the required rotational transform onl!.) Defining 
the function f(fj,i,Aip Py i) = - H-x^out - h), where 
ii nn and i> ou t are the transforms determined from solving 
Eq. (45) on the inner and outer interface for the mag- 
netic field parametrized by (/zj, A^> Pi j), and u-i and ti 
are prescribed input values, we employ a simple Newton 
method to set /(/i;, Aip Pt i) = 0. Typically, if a reason- 
able guess is provided, this converges in one or two itera- 
tions. Matrix perturbation methods are used to compute 
the derivatives: the infinitesimal variation, dA, resulting 
from an infinitesimal variation in Aipi or is given by 
A • d\ = dd — dA ■ A, and dA depends on A and dG 
via the chain rule; and similarly for dd. This search is 
computationally efficient: the integral metric elements, 
Eq. (37), do not need to be recomputed if the geometry 
does not change; and the matrix A is very sparse, and 
sparse linear algorithms are employed. 



pi 



1 0.05950 0.94168 (5 + 7 6)/(6 + 7 7) = 0.848898. 

2 0.35098 0.63872 (1 + 7 2)/(2 + 7 3) = 0.618034. 

3 0.64902 0.25740 (1 + 7 l)/(2 + 7 3) = 0.381966. 

4 1.00000 0.04106 (1 + 7 l)/(9 + 7 10) = 0.103971. 



In the MRXMHD model, the poloidal flux and helic- 
ity are assumed given. In the stepped-pressure model, 
the interface rotational transforms are constrained. (In 
both cases, the toroidal flux is constrained.) In either 
case, given the interface geometry, the Beltrami field in 
each volume that satisfies these constraints is unique - 
except for the possibility of bifurcations, which we do not 
consider in this article. So, with this implicit, the depen- 
dence of each Fi on the degrees of freedom is reduced to 
Fi=F t [x]. 

The task of constructing global equilibrium solutions 
is now the standard mathematical problem of finding ex- 
trema of the global energy functional, F = F[x\. Before 
describing the two basic approaches we have adopted for 
this, namely a preconditioned conjugate gradient algo- 
rithm for minimizing the global energy functional and a 
Newton-style algorithm for finding a zero of the multi- 
dimensional force-balance vector, we first present a con- 
vergence study illustrating that the Beltrami field in each 
Vi may be constructed to arbitrary accuracy. 



C. illustration of Beltrami fields 

For illustration, we show the Beltrami fields consis- 
tent with a multi-region equilibrium. The equilibrium is 
defined by the toroidal flux enclosed by each interface, 
the pressure in each Vi, and the interface transforms, 
and these are all given in Table. (I). The outer bound- 
ary is given by R = 1 + r cos(#) and Z = r sm(ff), with 
r '= 0.3 + 5cos{26 - () + <5cos(36> - Q and S = 1CT 3 . 
This choice of perturbation induces 'primary' islands at 
the * = 1/2 and t = 1/3 rational surfaces. The inte- 
rior boundaries are consistent with force balance. The 
interface cross sections and Poincare plots of the Bel- 
trami fields are shown in Fig. 1. Because of toroidal and 
poloidal coupling, magnetic islands (and irregular field 
lines) will form at all rational surfaces within the rota- 
tional transform range. 

With finite resolution, the equation V x B = /iB 
can of course only approximately be solved. However, 
given a smooth boundary, the solution to the Beltrami 
equation is well posed, and so the numerical error can 
be made arbitrarily small. Assuming that the Fourier 
resolution is sufficient to ensure the numerical error 
results from the finite-ness of the radial discretization, 
an n-th order approximation to the vector potential 
yields an error 0(h n+1 ), where h is the radial sub-grid 
size. The Fourier harmonics of the contravariant com- 



ponents of B are (^J~gB s 
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FIG. 1: Poincare plot showing the Beltrami fields in multiple, nested volumes and the ideal interfaces (thick lines) for the 
perturbed equilibrium on the cross sections (A) £ = 0; (B) £ = 7r/2; and (C) £ = tt. 
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FIG. 2: Scaling of components of error, 5j = j — ^iB, with respect to radial resolution. The diamonds are for the n = 3 (cubic) 
basis functions, the triangles are for the n = 5 (quintic) basis functions. The solid lines have gradient -3,-2 and -2, and the 
dotted lines have gradient -5,-4 and -4. 



(V9B 6 ) min = A' Ctmtn , and {^) m ^~A' 9m ^ 
where the prime denotes radial derivative. Radial 
derivatives reduce the order of the error, and B = V x A 
is generally an order less accurate that A itself; with 
the exception of B s , which remains accurate to 0(h n+1 ) 
as no radial derivatives are involved. Before computing 
V x B, the 'contravariant' components must be 'low- 
ered', B a = 9apB° , and the error in B s , Bg, B^ are 
each 0(h n ). The Fourier harmonics of the contravariant 
components of j = V x B arc computed similarly, and 
the error in j s , j 9 and f are <D(h n ), ©(/i™" 1 ) and 
0(/i" -1 ), respectively. The components of the error, 
j — /jB, arc quantified by 



/N_ 



AIN 



1/2 



for a = s, 6, and Q. These quantities are shown as a 
function of radial sub-grid resolution in Fig. 2, for the 
field in V3 of the equilibrium shown in Fig. 1. The 
expected error scalings, |<5j s | ~ 0(h n ), \Sj e \ ~ 0(ft™ -1 ) 
and \6j^\ ~ C(/i n_1 ), are confirmed for both the cubic, 
?i = 3, and quintic, n = 5, finite-element representations. 



Note that at no point does the algorithm depend of re- 
solving the fractal structure of phase space. The vector 
potential in each V; is a smooth function, both as a func- 
tion of position within a given volume, and as a function 
of interface geometry. 

Before proceeding to the task of piecing together mul- 
tiple, nested Beltrami fields to obtain global, non-trivial 
equilibria, we must tie down a 'loose-end' regarding the 
Fourier representation of the interfaces. 



D. spectral condensation 

To construct global equilibria, the Rij and Zij de- 
scribing the interface geometry will be varied to extrem- 
ize the energy functional and/or satisfy force balance. 
Tangential geometric variations merely change the angu- 
lar parametrization of the interfaces and do not change 
the interface geometry, and so do not affect the energy 
functional, but do alter the {Rij, Zij}. This freedom 
may be exploited to obtain a preferred angle parametriza- 
tion. 

A numerically insightful choice [7, 8] is to choose the 
angle that minimizes the 'spectral width', and so obtain 
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the most accurate representation of the interface geome- 
try for a given Fourier resolution. We define the spectral 
width as 



5 EC 



[Ri + zi 



(46) 



where p and q arc arbitrary integers required as input. 

The toroidal angle has already been constrained, 
( = — (f>, and the geometry of the interfaces is to 
be constrained by force balance, so we are left 
to minimize the spectral width with respect to 
poloidal variations, 5R = deRSu and SZ = dgZ 5u. To 
preserve stellarator symmetry we restrict attention 
to odd functions, Su = Uj sin(mj6 — rijQ. The 
variations in the Fourier harmonics, Rj and Zj, 
§ § RgSu cos(rrijO — UjQdOdC, and 



are given by 5Rj 



5Z j 
and Zg 
is 



ZgSusin(mj9 — UjC^dOdQ, where Rg = dgR 
dgZ . The first variation in the spectral width 



dOdC, (RgX + Z e Y) Su, (47) 



where X = X]j( m j + n j)Rj cos(m J -6 | — rijQ and 



Y = Ej( to j + n j)Zj sin(m J -0 - ,^ 
is decreased along 5u = —7, where I 
is extremized when 1 = 0. 



HjQ . The spectral width 
R e X + Z 8 Y, and 



E. illustrations of global equilibria 

We have implemented two numerical meth- 
ods for finding global equilibria. The first is 
a minimization algorithm: we seek minima of 
-H x ] = J2i Jybi/(7 - 1) + B 2 /2]dv using a precon- 
ditioned, conjugate gradient algorithm, where in each 
volume piV 1 = ai is constant, and the field satisfies 
V x B = /i;B. The gradient of F with respect to the 
Rij and Zij, with the fluxes and helicity constrained, 
can be derived from Eq. (12) by recalling that the 
displacement is £ = 8R R + 6Z k and the surface element 
is ds = -RZ R + {Z g R c - RgZ c )cj) + RRgz. The 
derivatives of F are 



dF 
dF 



dZ, 



1,3 



= - ([]p + B 2 /2]]iRZ ). - (Reli)j, (48) 
= + {[\p + B 2 /2]] l RR 6 ). - (Zgh) 3 , (49) 



where the (Rgli)j and (Zgli)j terms are included to re- 
duce the spectral width. 

The second approach for finding global equilibria is a 
globally-convergent, multi-dimensional Newton method. 
We seek a zero of the 'force-balance' vector, constructed 
as follows. Global force balance is satisfied when the to- 
tal pressure discontinuity across each interface is zero, 
[[p + £> 2 /2]] = 0. The Fourier representation of the inter- 
faces minimizes the spectral width when / is zero. Thus, 
we construct a 'constraint' vector, f(x), by collecting to- 
gether the harmonics [[p + Ft 2 /2}]i,j and Ii.j. 

Within the stellarator symmetric representation we 
have employed, R is an even function of (9, £), an d Z 



0.3 ■ 



0.2 



N 0.0 ■ 



-0.1 - 



-0.2 - 



-0.3 ■ 



0.7 




' .3 



FIG. 3: Comparison between the SPEC interfaces, with 
Nv = 6, and the corresponding VMEC surfaces (thick lines, 
upper half); and the SPEC radial sub-grid (lower half). 
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FIG. 4: Location of magnetic axis as computed by SPEC 
against resolution, Nv, of stepped-pressure approximation to 
smooth pressure profile. The dotted line is the location of 
magnetic axis for a high resolution VMEC calculation. 



is odd. Force-balance, [[p + £? 2 /2]], and the spectral min- 
imization condition, /, are similarly even and odd func- 
tions. So, after suitably truncating [[p + B 2 /2]} and / 
to match the truncated Fourier representation of the in- 
terface geometry, the number of constraints equals the 
number of degrees-of-freedom. 

Expanding about an arbitrary point, 
f (x + <5x) sa f (x) + Vf (x) • (5x, the Newton correc- 
tion required to find the equilibrium point, f = 0, is 
given by Sx = — Vf -1 • f. The Jacobian matrix, Vf, 
describes how the Beltrami fields change when the geom- 
etry is changed (more precisely, how B 2 on the interfaces 
changes when the I; change) and is computed in parallel 
using matrix perturbation methods as described above. 

In the following, we use this method for constructing 
global solutions. In all calculations presented, either ex- 
plicitly or implicitly, for any given Fourier and radial 
sub-grid resolution, the error in force balance, |/|, and 
the error in 'position', |<5x|, are less than 10~ 12 . 

Before presenting illustrations of non-axisymmetric 
global equilibria, we first present a comparison of 
stepped-pressure equilibria to axisymmctric MHD equi- 
libria with smooth profiles, the latter computed by 
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smooth pressure profile. 



VMEC [84]. SPEC is intended for computing equilib- 
ria with partially chaotic fields and only admits stepped- 
pressure profiles. In contrast, VMEC globally constrains 
the field to be integrable, and so admits (and assumes) 
smooth profiles. As it is the profiles that define an equi- 
librium, VMEC and SPEC will differ. To obtain agree- 
ment, it is required to perform a convergence study in the 
pressure profile: as the number of steps in the stepped 
profile increases, the stepped profile will better approxi- 
mate a smooth profile. 

We consider an equilibrium with boundary described 
by R(6,Q = 1.0 + 0.3 cosfl and Z(6, Q = 0.3 sin0. 
For the (smooth) pressure profile we take 
p(ip) — Po(l ~ 2?/> + -0 2 ), where ip is the normalized 
toroidal flux and po is a scaling factor chosen to give a 
Shafranov shift about one-third the minor radius. For 
the (smooth) transform profile we take * = *o — *iVa 
where t = (8 + 97)/ (9 + IO7) and t x = t Q - (1 + I7)/ 
(9 + IO7). We use high Fourier resolution, M = 16, and 
high radial sub-grid resolution; this will ensure that any 
discrepancy results from finite Ny, and so the error will 
decrease as Ny is increased. 

As input to SPEC, we must 'discrctize' the pressure 
and transform profiles. Earlier we recalled that invari- 
ant surfaces of Hamiltonian systems with noble rotational 
transform are the most likely to survive chaos-inducing 
perturbations, and so the interface transforms should be 
strongly irrational; however, in the axisymmctric and 
therefore integrable case, the interface transforms may be 
chosen arbitrarily. For convergence studies, it is prefer- 
able to have interfaces regularly spaced in radius and 
we choose the interface rotational transforms */ = t(ipi), 
where y^; = l/Ny. We discretize the pressure profile in 
such a way as to piecewise conserve the integrated pres- 
sure, 

ripi ri>i 
pi I dip = I p(tp) dtp. (50) 
Jil>i-i J'ipi-i 

For illustration, the SPEC interfaces and the corre- 
sponding surfaces (identified by toroidal flux) of a high 
radial resolution VMEC equilibrium are shown for a low 
Ny case in Fig. 3. Despite the fact that a smooth pro- 
file was supplied to VMEC but a stepped-prcssure profile 
was supplied to SPEC, the agreement is reasonable. This 
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FIG. 6: Difference between finite M, N approximation to 
interface geometry, and a high-resolution reference approxi- 
mation (with M = 13 and N = 8), plotted against Fourier 
resolution. 



may be expected: the Shafranov shift depends [95] pri- 
marily on the integral of the pressure profile and not, to 
lowest order, on the precise details of the pressure pro- 
file itself. To quantify the discrepancy, we compare the 
location of the magnetic axes. As shown in Fig. 4, as 
Ny is increased the agreement improves. The stepped 
pressure profile approximation and the smooth pressure 
profile are shown in Fig. 5. 

Indeed, it may be shown that the infinite interface limit 
of MRXMHD (i.e. assuming continuously nested mag- 
netic flux surfaces) reduces to Vp = j x B and j • n = 0, 
where n is normal to the flux surfaces. These are the 
equations that define ideal MHD equilibria. The details 
of this derivation will be presented in a forthcoming pa- 
per. 

To illustrate a non-axisymmctric global equilibrium, 
we return to the 'perturbed' equilibrium shown in Fig. 1. 
Assuming that the radial sub-grid resolution is suffi- 
ciently high so that the error is dominated by finite 
Fourier resolution, we now confirm that the error in the 
interface geometry decreases as M and N increase. 

Let the exact solution for the Z-th interface be de- 
scribed by R(8,Q and Z(9,Q- The error between 
this and an approximation, Rm.n(<x,C) an d Zm,n{oi,(,), 
with a potentially different poloidal angle, a, on the 
£0 plane is quantified by A = J D(9) dl, where dl 2 = 
dx 2 + dy 2 is the arc-length of the curve x{9) = R(9, Co) 
and y(9) = Z(9,(o), and D is the distance between 
this reference curve and the curve described by 
xm,n{cx) = Rm,n{ol, Co) and j/M,Jv(a) = Zm,n(cx, Co), i-e. 
D 2 = [x(6)-x u ,N(a)} 2 + [y(e)-y M ,N(a)} 2 . The 
comparison is made at the same polar angle, so that 
y(9)/[x{9) - x ] = yM.N{a)/[xM,N{a) - x ], where x is 
a reference point (e.g. the magnetic axis). 

As the exact solution is not known apriori, we take as 
the reference configuration the highest resolution approx- 
imation available. The error in the interface geometry, for 
each of the internal interfaces, is shown as a function of 
(M,N) in Fig. 6, where we see that the error decreases 
as the Fourier resolution is increased. 

All properties of the equilibrium are defined by the 
interface geometries: if the interface geometry has con- 
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verged, then so too have the Beltrami fields, and the lo- 
cation and size of the magnetic islands and chaotic seas 
- well, it is difficult to resolve the infinitely complicated 
structure of phase space, but B itself and all integral 
properties are converged. 

The good convergence properties of the interface ge- 
ometry with Fourier resolution is because (i) the inter- 
faces are chosen to have the most noble transform, and 
so are 'as far away as possible', so to speak, from the 
lowest order islands and associated chaos, and so are 
the smoothest surfaces (this is in contrast to flux sur- 
faces adjacent to a separatrix, or flux surfaces that are 
nearly critical); and (ii) that the Fourier representation 
exploits a preferred angle parametrization that minimizes 
the spectral width. More importantly, perhaps, is that 
convergence is obtained because there is a well defined, 
exact solution, and that the numerical discretization is 
capable of resolving all the structure of the solution. 

As a final illustration, we present a stepped-pressure 
equilibrium consistent with the boundary and profiles ob- 
tained via a 3D STELLOPT reconstruction [96] of an up- 
down symmetric DIIID experimental shot with applied 
resonant magnetic perturbation (RMP) fields. The re- 
construction process seeks to infer the experimental con- 
figuration by adjusting the MHD equilibrium (presently, 
STELLOPT is built around VMEC) by varying the the 
plasma boundary and the pressure and current profiles, 
so that derived quantities (such as Thomson scattering, 
motional Stark effect polarimetry, and magnetic diagnos- 
tics) match the experimental measurements. Because of 
the applied error fields, and the plasma response to these 
error fields, the reconstructed boundary is slightly, but 
significantly, perturbed from axisymmctry. The pressure 
and g-profiles, where q is the safety- factor q = 1/t, de- 
rived from the reconstruction are shown in Fig. 7. It 
is interesting to observe that the reconstructed pressure 
profile appears quite flat across the lowest order ratio- 
nal surfaces. Furthermore, the locations of locally maxi- 
mum pressure gradient appear to coincide with strongly- 
irrational surfaces. 

We compute the stepped-pressure equilibrium using 
the reconstructed boundary, a stepped, Ny = 32 approx- 
imation (Fig. 7) to the reconstructed pressure profile, and 
the reconstructed g-profile. The rotational transforms of 
the interfaces are chosen by selecting the most noble ir- 
rationals that are within range. The Fourier resolution 
is M = 10 and N = 6, and the total radial sub-grid 
resolution is 279. A Poincare plot is shown in Fig. 8; 
most visible is a q = 2 island at where VMEC has the 
q = 2 rational surface. In this 'fixed-boundary' calcula- 
tion, the boundary is constrained to remain a fixed, good 
flux surface. To determine to what extent the RMP fields 
and the plasma response ergodize the field in the vicin- 
ity of the plasma edge, a 'free-boundary' calculation is 
required. This is left for future work. 



V. COMMENTS, FUTURE WORK 

The fact that stepped-pressure equilibria can be de- 
rived as minima of an energy functional is a great con- 




FIG. 7: Pressure profile (smooth) from a DIIID reconstruc- 
tion using STELLOPT and stepped-pressure approximation. 
Also shown is the inverse rotational transform = safety factor. 
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FIG. 8: Poincare plot of a DIIID equilibrium with perturbed 
boundary, calculated using SPEC. 



venience numerically, as this allows employment of mini- 
mization methods to construct the Beltrami fields and 
interface geometries that satisfy force balance. Fur- 
thermore, the MRXMHD energy functional provides a 
self-consistent approach for determining the stability of 
partially chaotic equilibria [75-77]. In future work, we 
hope to explore whether the suppression of edge local- 
ized modes by resonant magnetic perturbations [16], that 
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result in a stochastic plasma edge, may be understood 
through an MRXMHD stability analysis. 

In the above text we have distinguished MRXMHD 
equilibria, for which the algorithm allows fj, to vary to 
force the helicity constraint, from stepped-pressure equi- 
libria, for which the algorithm allows both \i and the 
poloidal flux to vary to force the transform constraint and 
for which the self-consistent helicity is computed apostc- 
ori. The distinction is, however, superficial. The profiles 
could also be specified by keeping \i and the poloidal flux 
as fixed input parameters. That the profiles can be pre- 
scribed and constrained in a variety of ways illustrates 
the flexibility of the theoretical and numerical method. 
(Note that because ideal MHD is not applied globally, 
there is no explicit relationship between the toroidal and 
poloidal fluxes and the rotational-transform profile, for 
example.) To decide how the profiles should be described, 
one must depart from equilibrium theory and develop 
a self-consistent model of transport, which may suggest 
how, for example, the poloidal flux, the parallel current 
and the helicity should vary in time to preserve the trans- 
form constraint. 

Given that the MRXMHD equilibria do not have 1/x 
and (5-function currents at the rational surfaces, and that, 
in a pressure transport model that includes a small per- 
pendicular diffusion, maximum pressure gradients will 
appear at the most irrational locations [32, 97] - and 
so too will, perhaps, the pressure driven currents - then 
it would seem that MRXMHD equilibria are smoothly 
connected (e.g. as the perpendicular-diffusion coefficient 
k± —> 0) to nearly-ideal steady-state equilibria, and so it 
may be advantageous to initialize resistive MHD initial 
value codes such as NIMROD [51] or M3DC1 [50] with 
SPEC equilibria. Local flattening of the pressure gradi- 
ent across the resonances can provide increased stability 
and can allow access to higher plasma beta [98, 99]. 

The equilibrium model considered in this article docs 
not include plasma flow, which can impact both the 
equilibrium and island healing phenomena [100, 101]. 
Extended MHD modeling of plasmas [98, 99, 102-104] 
is crucial for understanding some key experimental ob- 
servations. Perhaps, by extending a stepped-pressure 
equilibrium, an equilibrium with a continuous-pressure 
profile could be constructed by replacing the pressure 
jump interfaces with finite-width ideal-MHD layers [76], 
each topologically constrained to avoid resonances and 
to avoid ideal instabilities [76], and, perhaps, capable 
of supporting extended MHD behavior. A variational 
principle based on minimizing the generalized enstrophy 
[105], F = /|Vx (V + A)\ 2 dx, may better describe self- 



organization in two-fluid plasmas. 

In addition to the incremental code improvements that 
are inevitably required, the following in particular will be 
pursued. In this paper it was assumed that the plasma 
boundary, dV, is a given, fixed toroidal surface. More 
generally, the magnetic field is part generated by external 
currents, and the /ree-boundary problem could be solved 
with a little extra work, where dV is to be determined 
as part of the solution. Also, most modern tokamaks are 
not up-down symmetric, so it will be required to relax 
the stellar ator symmetry constraint. 

It is worth exploring more efficient numerical methods 
for computing the Beltrami fields. In a closed domain 
V in R 3 , in general multiply connected, the solutions of 
VxB = /iB can be represented [38, 39] by 

B = (Vx+/*)/ G(r,r')B'xn'dS', (51) 

JdV 

where n is the outward unit normal on dV and G(r, r') 
satisfies (V 2 + /i 2 )G(r, r') = — S(r — r'), with <$(•) being 
the 3D Dirac 5 function. 

Given the geometry of an interface, the maximum pres- 
sure jump that an interface can support can be quickly 
determined by an analysis of the pressure-jump Hamil- 
tonian: the pressure discontinuity, 2(p + — p~), is in- 
creased until the appropriate irrational surface of the 
pressure-jump Hamiltonian is critical. Looking beyond 
our present task of constructing an equilibrium consis- 
tent with a given pressure that is not changed, this gives 
an efficient method for distributing the pressure so that 
the most robust interfaces support the most pressure. 
As the pressure across any interface is altered, there will 
be a global response that requires re-computation of the 
equilibrium. 

We dedicate this article to Paul Garabedian, whom we 
consider to be a pioneer in the field of 3D MHD cal- 
culations, and who also endorsed and employed weak 
solutions [106]. We are also indebted to Steve Hir- 
shman, who has provided VMEC. One of us (SRH) 
acknowledges stimulating discussions with Neil Pom- 
phrey, with Don Monticello regarding the construction 
of straight-ficld-linc coordinates, and with Allen Boozer 
regarding the choice of gauge [107] for the vector po- 
tential, and RLD acknowledges discussion with Robert 
MacKay. The authors gratefully acknowledge support of 
the U.S. Department of Energy and the Australian Re- 
search Council, through grants DP0452728, FT0991899 
and DP110102881. 
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